Categoría: Sol_1Eva 2011-2020

  • s1Eva_IIT2018_T2 Distancia mínima a un punto

    Ejercicio: 1Eva_IIT2018_T2 Distancia mínima a un punto

    Literal a

    Se requiere analizar la distancias entre una trayectoria y el punto = [1,1]

    Al analizar las distancias de ex y el punto [1,1] se trazan lineas paralelas a los ejes desde el punto [1,1], por lo que se determina que el intervalo de x = [a,b] para distancias se encuentra en:

    a > 0, a = 0.1
    b < 1, b = 0.7

    El ejercicio usa la fórmula de distancia entre dos puntos:

    d = \sqrt{(x_2-x_1)^2+(y_2- y_1)^2}

    en los cuales:

    [x1,y1] = [1,1]
    [x2,y2] = [x, ex]

    que al sustituir en la fórmula se convierte en:

    d = \sqrt{(x-1)^2+(e^x- 1)^2}

    que es lo requerido en el literal a


    Literal b

    Para usar un método de búsqueda de raíces, se requiere encontrar el valor cuando f(x) = d' = 0.

    Un método como el de Newton Raphson requiere también f'(x) = d''

    f(x) = \frac{x + (e^x - 1)e^x - 1}{\sqrt{(x - 1)^2 + (e^x - 1)^2}} f'(x)= \frac{(e^x - 1)e^x + e^{2x} + 1 - \frac{(x + (e^x - 1)e^x - 1)^2}{(x - 1)^2 + (e^x - 1)^2}} {\sqrt{(x - 1)^2 + (e^x - 1)^2}}

    expresiones obtenidas usando Sympy

    f(x) :
    (x + (exp(x) - 1)*exp(x) - 1)/sqrt((x - 1)**2 + (exp(x) - 1)**2)
    f'(x) :
    ((exp(x) - 1)*exp(x) + exp(2*x) + 1 - (x + (exp(x) - 1)*exp(x) - 1)**2/((x - 1)**2 + (exp(x) - 1)**2))/sqrt((x - 1)**2 + (exp(x) - 1)**2)
    
    f(x) :
           / x    \  x        
       x + \e  - 1/*e  - 1    
    --------------------------
        ______________________
       /                    2 
      /         2   / x    \  
    \/   (x - 1)  + \e  - 1/  
    f'(x) :
                                                  2
                             /    / x    \  x    \ 
    / x    \  x    2*x       \x + \e  - 1/*e  - 1/ 
    \e  - 1/*e  + e    + 1 - ----------------------
                                                 2 
                                     2   / x    \  
                              (x - 1)  + \e  - 1/  
    -----------------------------------------------
                   ______________________          
                  /                    2           
                 /         2   / x    \            
               \/   (x - 1)  + \e  - 1/            
    
    
    

    lo que permite observar la raíz de f(x) en una gráfica:
    distancia mínima f(x)
    con las siguientes instrucciones:

    # Eva_IIT2018_T2 Distancia mínima a un punto
    import numpy as np
    import matplotlib.pyplot as plt
    import sympy as sym
    
    # INGRESO
    x = sym.Symbol('x')
    fx = sym.sqrt((x-1)**2+(sym.exp(x) -1)**2)
    
    a = 0
    b = 1
    muestras = 21
    
    # PROCEDIMIENTO
    dfx = sym.diff(fx,x,1)
    d2fx = sym.diff(fx,x,2)
    
    f = sym.lambdify(x,dfx)
    xi = np.linspace(a,b,muestras)
    fi = f(xi)
    
    
    # SALIDA
    print('f(x) :')
    print(dfx)
    print("f'(x) :")
    print(d2fx)
    print()
    print('f(x) :')
    sym.pprint(dfx)
    print("f'(x) :")
    sym.pprint(d2fx)
    
    # GRAFICA
    plt.plot(xi,fi, label='f(x)')
    plt.axhline(0)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.legend()
    plt.grid()
    plt.show()
    

    Usando el método de la bisección para el intervalo dado, se tiene:

    f(x) = \frac{x + (e^x - 1)e^x - 1}{\sqrt{(x - 1)^2 + (e^x - 1)^2}}

    itera = 0 , a = 0, b=1

    c= \frac{0+1}{2} = 0.5 f(0) = \frac{0 + (e^0 - 1)e^0 - 1}{\sqrt{(0 - 1)^2 + (e^0 - 1)^2}} = -1 f(1) = \frac{1 + (e^1 - 1)e^1 - 1}{\sqrt{(1 - 1)^2 + (e^1 - 1)^2}} 2.7183 f(0.5) = \frac{(0.5) + (e^(0.5) - 1)e^(0.5) - 1}{\sqrt{((0.5) - 1)^2 + (e^(0.5) - 1)^2}} = 0.6954

    cambio de signo a la izquierda,

    a= 0, b=c=0.5

    tramo = |0.5-0| = 0.5

    itera = 1

    c= \frac{0+0.5}{2} = 0.25 f(0.25) = \frac{(0.25) + (e^(0.25) - 1)e^(0.25) - 1}{\sqrt{((0.25) - 1)^2 + (e^(0.25) - 1)^2}} = -0.4804

    cambio de signo a la derecha,

    a=c= 0.25, b=0.5

    itera = 2

    c= \frac{0.25+0.5}{2} = 0.375 f(0.375) = \frac{(0.375) + (e^(0.375) - 1)e^(0.375) - 1}{\sqrt{((0.375) - 1)^2 + (e^(0.375) - 1)^2}} = 0.0479

    cambio de signo a la izquierda,

    a= 0.25, b=c=0.375

    se continúan las iteraciones con el algoritmo, para encontrar la raíz en 0.364:

    método de Bisección
    i ['a', 'c', 'b'] ['f(a)', 'f(c)', 'f(b)']
       tramo
    0 [0, 0.5, 1] [-1.      0.6954  2.7183]
       0.5
    1 [0, 0.25, 0.5] [-1.     -0.4804  0.6954]
       0.25
    2 [0.25, 0.375, 0.5] [-0.4804  0.0479  0.6954]
       0.125
    3 [0.25, 0.3125, 0.375] [-0.4804 -0.2388  0.0479]
       0.0625
    4 [0.3125, 0.34375, 0.375] [-0.2388 -0.1004  0.0479]
       0.03125
    5 [0.34375, 0.359375, 0.375] [-0.1004 -0.0274  0.0479]
       0.015625
    6 [0.359375, 0.3671875, 0.375] [-0.0274  0.01    0.0479]
       0.0078125
    7 [0.359375, 0.36328125, 0.3671875] [-0.0274 -0.0088  0.01  ]
       0.00390625
    8 [0.36328125, 0.365234375, 0.3671875] [-0.0088  0.0006  0.01  ]
       0.001953125
    9 [0.36328125, 0.3642578125, 0.365234375] [-0.0088 -0.0041  0.0006]
       0.0009765625
    raíz en:  0.3642578125
    

    Al algoritmo anterior se complementa con las instrucciones de la función para la bisección.

    # Eva_IIT2018_T2 Distancia mínima a un punto
    import numpy as np
    import matplotlib.pyplot as plt
    import sympy as sym
    
    # INGRESO
    x = sym.Symbol('x')
    fx = sym.sqrt((x-1)**2+(sym.exp(x) -1)**2)
    
    a = 0
    b = 1
    muestras = 21
    
    # PROCEDIMIENTO
    dfx = sym.diff(fx,x,1)
    d2fx = sym.diff(fx,x,2)
    
    f = sym.lambdify(x,dfx)
    xi = np.linspace(a,b,muestras)
    fi = f(xi)
    
    
    # SALIDA
    print('f(x) :')
    print(dfx)
    print("f'(x) :")
    print(d2fx)
    print()
    print('f(x) :')
    sym.pprint(dfx)
    print("f'(x) :")
    sym.pprint(d2fx)
    
    # GRAFICA
    plt.plot(xi,fi, label='f(x)')
    plt.axhline(0)
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.legend()
    plt.grid()
    plt.show()
    
    # Algoritmo de Bisección
    # [a,b] se escogen de la gráfica de la función
    # error = tolera
    import numpy as np
    
    def biseccion(fx,a,b,tolera,iteramax = 20, vertabla=False, precision=4):
        '''
        Algoritmo de Bisección
        Los valores de [a,b] son seleccionados
        desde la gráfica de la función
        error = tolera
        '''
        fa = fx(a)
        fb = fx(b)
        tramo = np.abs(b-a)
        itera = 0
        cambia = np.sign(fa)*np.sign(fb)
        if cambia<0: # existe cambio de signo f(a) vs f(b)
            if vertabla==True:
                print('método de Bisección')
                print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
                print('  ','tramo')
                np.set_printoptions(precision)
                
            while (tramo>=tolera and itera<=iteramax):
                c = (a+b)/2
                fc = fx(c)
                cambia = np.sign(fa)*np.sign(fc)
                if vertabla==True:
                    print(itera,[a,c,b],np.array([fa,fc,fb]))
                if (cambia<0):
                    b = c
                    fb = fc
                else:
                    a = c
                    fa = fc
                tramo = np.abs(b-a)
                if vertabla==True:
                    print('  ',tramo)
                itera = itera + 1
            respuesta = c
            # Valida respuesta
            if (itera>=iteramax):
                respuesta = np.nan
    
        else: 
            print(' No existe cambio de signo entre f(a) y f(b)')
            print(' f(a) =',fa,',  f(b) =',fb) 
            respuesta=np.nan
        return(respuesta)
    
    # INGRESO
    tolera = 0.001
    
    # PROCEDIMIENTO
    respuesta = biseccion(f,a,b,tolera,vertabla=True)
    # SALIDA
    print('raíz en: ', respuesta)
    

    .

  • s1Eva_IT2018_T1 Tanque esférico canchas deportivas

    Ejercicio: 1Eva_IT2018_T1 Tanque esférico canchas deportivas

    a) Para seleccionar el rango para h=[a,b], se observa que el tanque puede estar vacío, a=0 o lleno al máximo, b=2R = 2(3)=6, con lo que se obtiene:

    h =[0.0, 6.0]

    conociendo la proporción con el valor máximo, se tiene un valor inicial para h0 para las iteraciones.

    Vmax = \frac{\pi}{3} (2R)^2 (3R-2R) Vmax = \frac{4\pi }{3}R^{3}= 113.09 h_0 = (6)*30/113.09 = 1.59

    b) Usar Newton-Raphson con tolerancia 1e-6

    f(h) = \frac{\pi }{3}h^2 (3(3)-h)-30 f(h) = \frac{\pi }{3}(9h^2 -h^3-28.647) f'(h) = \frac{\pi }{3}(18h-3h^2)

    el punto siguiente de iteración es:

    h_{i+1} = h_{i} -\frac{f(h)}{f'(h)} = h_{i}-\frac{ \frac{\pi }{3}(9h^2 -h^3-28.647)}{ \frac{\pi }{3}(18h-3h^2)} h_{i+1} = h_{i} -\frac{(9h^2 -h^3-28.647)}{(18h-3h^2)}

    con lo que se realiza la tabla de iteraciones:

    hi hi+1 error orden
    1.590 2.061 0.47 10-1
    2.061 2.027 -0.034 10-2
    2.027 2.02686 -0.00014 10-4
    2.02686 2.0268689 -2.32E-09 10-9

    En valor práctico 2.028 m usando flexómetro, a menos que use medidor laser con precisión 10-6 usará más dígitos con un tanque de 6 metros de altura ganará una precisión de una gota de agua para usar en duchas o regar el césped .

    c) El orden de convergencia del error observando las tres primeras iteraciones es cuadrático

    Tarea: Realizar los cálculos con Python, luego aplique otro método. Añada sus observaciones y conclusiones.

  • s1Eva_IT2018_T3 Temperatura en nodos de placa

    Ejercicio: 1Eva_IT2018_T3 Temperatura en nodos de placa

    a) Plantear el sistema de ecuaciones. Usando el promedio para cada nodo interior: Placa Temp 03

    a=\frac{50+c+100+b}{4} b=\frac{a+30+50+d}{4} c=\frac{a+60+100+d}{4} d=\frac{b+60+c+30}{4}

    que reordenando se convierte en:

    4a=150+c+b 4b=a+80+d 4c=a+160+d 4d=b+c+90

    simplificando:

    4a-b-c= 150 a-4b+d = -80 a-4c+d = -160 b+c-4d = -90

    que a forma matricial se convierte en:

    A = [[ 4, -1, -1, 0.0],
         [ 1, -4,  0, 1.0],
         [ 1,  0, -4, 1.0],
         [ 0,  1,  1,-4.0]]
    B = [[ 150.0],
         [ -80.0],
         [-160.0],
         [ -90.0]]
    

    Observación: la matriz A ya es diagonal dominante, no requiere pivotear por filas.  Se aumentó el punto decimal a los valores de la matriz A y el vector B  para que sean considerados como números reales.

    El número de condición es: np.linalg.cond(A) = 3.0

    que es cercano a 1 en un orden de magnitud, por lo que la solución matricial es "estable" y los cambios en los coeficientes afectan proporcionalmente a los resultados. Se puede aplicar métodos iterativos sin mayores inconvenientes.

    b y c) método de Jacobi para sistemas de ecuaciones, con vector inicial

     
    X(0) = [[60.0],
            [40], 
            [70],
            [50]] 
    

    reemplazando los valores iniciales en cada ecuación sin cambios.

    iteración 1
    a=\frac{50+70+100+40}{4} = 65

    b=\frac{60+30+50+50}{4} = 47.5 c=\frac{60+60+100+50}{4} = 67.5 d=\frac{40+60+70+30}{4} = 50
    X(1) = [[65],
            [47.5],
            [67.5],
            [50]]
    
    vector de error = 
         [|65-60|,
          |47.5-40|,
          |67.5-70|,
          |50-50|]
      = [|5|,
         |7.5|,
         |-2.5|,
         |0|]
    errormax = 7.5
    

    iteración 2
    a=\frac{50+67.5+100+47.5}{4} = 66.25

    b=\frac{65+30+50+50}{4} = 48.75 c=\frac{65+60+100+50}{4} = 68.75 d=\frac{47.5+60+67.7+30}{4} = 51.3
    X(2) = [[66.25],
            [48.75],
            [68.75],
            [51.3]]
    
    vector de error = 
           [|66.25-65|,
            |48.75-47.5|,
            |68.75-67.5|,
               |51.3-50|] 
      = [|1.25|,
         |1.25|,
         |1.25|,
         |1.3|]
    errormax = 1.3
    

    iteración 3
    a=\frac{50+68.75+100+48.75}{4} = 66.875

    b=\frac{66.25+30+50+51.3}{4} = 49.38 c=\frac{66.25+60+100+51.3}{4} = 69.3875 d=\frac{48.75+60+68.75+30}{4} = 51.875
    X(2) = [[66.875],
            [49.38],
            [69.3875],
            [51.875]]
    
    vector de error = 
          [|66.875-66.25|,
           |49.38-48.75|,
           |69.3875-68.75|,
           |51.875-51.3|]
        = [|0.655|,
           |0,63|,
           |0.6375|,
            |0.575|]
    errormax = 0.655
    con error relativo de:
    100*0.655/66.875 = 0.97%

    siguiendo las iteraciones se debería llegar a:

    >>> np.linalg.solve(A,B)
    array([[ 67.5],
           [ 50. ],
           [ 70. ],
           [ 52.5]])
    
  • s1Eva_IT2018_T4 El gol imposible

    Ejercicio: 1Eva_IT2018_T4 El gol imposible

    Tabla de datos:

    ti = [0.00, 0.15, 0.30, 0.45, 0.60, 0.75, 0.90, 1.05, 1.20]
    xi = [0.00, 0.50, 1.00, 1.50, 1.80, 2.00, 1.90, 1.10, 0.30]
    yi = [0.00, 4.44, 8.88,13.31,17.75,22.19,26.63,31.06,35.50]
    zi = [0.00, 0.81, 1.40, 1.77, 1.91, 1.84, 1.55, 1.03, 0.30]
    

    Observe que, un gol simplificado con física básica es un tiro parabólico cuya trayectoria se compone de movimientos en los ejes, Y y Z. Sin embargo, lo "imposible" del gol mostrado es añadir el movimiento en X. (Referencia de la imagen en el enunciado)

    El movimiento de "profundidad" o dirección hacia el arco y(t) es semejante a un polinomio de primer grado, y el movimiento de "altura" z(t) es un polinomio de segundo grado. El movimiento "imposible" en el eje X, podría ser descrito por un polinomio de segundo o mayor grado.

    a) Encontrar t para altura máxima, que se encuentra al igualar la derivada dz/dt a cero. Para interpolar el polinomio z(t), de segundo grado, se puede usar tres puntos de los sugeridos: 0, 0.3 y 0.6, que además son equidistantes en t (facilita usar cualquier método de interpolación).

    Por ejemplo, con diferencias finitas avanzadas:

    t z(t) d1 d2 d3
    0.00 0.00 1.40 -0.89
    0.30 1.40 0.51
    0.60 1.91
    z(t) = 0 + 1.40\frac{(t-0)}{1!(0.3)} - 0.89 \frac{(t-0)(t-0.3)}{2!(0.3)^2} = 0 + 4.66 t - 4.94(t^2-0.3t) = 4.66 t + 1.48 t - 4.94 t^2 z(t) = 6.42 t - 4.94 t^2

    para encontrar el valor máximo se encuentra \frac{dz}{dt} = 0

    \frac{dz}{dt} = 6.42 - 2(4.94) t 6.42 - 2(4.94) t = 0 t = \frac{6.42}{2(4.94)} t = 0.649

    Observe que el resultado tiene sentido, pues según la tabla, el máximo debería estar entre 0.60 y 0.75

    b) El cruce por la "barrera", corresponde al desplazamiento del balón en el eje Y a 9 metros: y(t) = 9.
    El polinomio modelo de primer grado usa dos puntos para la interpolación, de los sugeridos se escogen dos, por ejemplo: 0.0 y 0.3.

    Usando diferencias finitas avanzadas :

    d1 = (8.88-0) = 8.88 y(t) = 0 + 8.88\frac{(t-0)}{1!(0.3)} y(t) = 29.6 t

    usando y(t) = 9

    29.6 t = 9
    t = 0.30
    z(0.30) = 1.40
    (de la tabla de datos)

    cuya respuesta con solo dos dígitos decimales es coherente, al estar cercano el valor a una distancia y=8.88 o aproximado a 9.
    La respuesta puede ser mejorada usando más digitos significativos en las operaciones.

    c)  La desviación máxima en eje X.
    Determine un polinomio para la trayectoria en el eje X y obtenga el máximo igualando la derivada del polinomio x(t) a cero.

    Por simplicidad, se usa un polinomio de segundo grado, alrededor de los valores máximos en el eje X

    t x(t) d1 d2 d3
    0.60 1.80 0.20 -0.30
    0.75 2.00 -0.10
    0.90 1.90
    x(t) = 1.80 + 0.20 \frac{(t-0.60)}{1!(0.15)} -0.30 \frac{(t-0.60)(t-0.75)}{2!(0.15)^2} x(t) = 1.80 + 1.33 (t-0.60) - 6.67(t-0.60)(t-0.75)

    como se busca el máximo, se usa la derivada \frac{dx}{dt} = 0

    \frac{dx}{dt} = 1.33 - 6.67(2t +(-0.60-0.75)) 1.33 - 13.34t + 9.00 = 0 -13.34t + 10.33 = 0

    t = 0.77

    x(0.77) = 1.80 + 1.33(0.77-0.60) - 6.67(0.77-0.60)(0.77-0.75) x(0.77) = 2.003

    lo que es coherente con la tabla para el eje x, pues el máximo registrado es 2, y el próximo valor es menor, la curva será decreciente.


    Desarrollo extra para observar y verificar resultados:

    Usando los puntos y un graficador 3D se puede observar la trayectoria:

    Gol Imposible 02

    Tarea: Realice el ejercicio usando los algoritmos en Python, muestre los polinomios obtenidos y compare.

    Nota: La gráfica 3D presentada usa interpolación de Lagrange con todos los puntos. Realice las observaciones y recomendaciones necesarias y presente su propuesta como tarea.

  • s1Eva_IIT2017_T1 Aproximar a polinomio usando puntos

    Ejercicio: 1Eva_IIT2017_T1 Aproximar a polinomio usando puntos

    Se dispone de tres puntos para la gráfica.

    x  f(x)
     0  1
     0.2  1.6
     0.4  2.0

    Si el polinomio de Taylor fuera de grado 0, sería una constante, que si se evalúa en x0 = 0 para eliminar los otros términos, se encuentra que sería igual a 1

    Como se pide el polinomio de grado 2, se tiene la forma:

    p(x) = a + bx + c x ^2 p(x) = 1 + bx + c x^2

    Se disponen de dos puntos adicionales que se pueden usar para determinar b y c:

    p(0.2) = 1 + 0.2 b + (0.2)^2 c = 1.6 p(0.4) = 1 + 0.4 b + (0.4)^2 c = 2.0

    simplificando:

    0.2 b + (0.2)^2 c = 1.6 - 1 = 0.6 0.4 b + (0.4)^2 c = 2.0 - 1 = 1

    multiplicando la primera ecuación por 2 y restando la segunda ecuación:

    0 - 0.08 c = 1.2-1 = 0.2 c = - 0.2/0.08 = -2.5

    sustituyendo el valor de c obtenido en la primera ecuación

    0.2 b + 0.04(-2.5) = 0.6 0.2 b = 0.6 - 0.04(-2.5) = 0.6 + 0.1 = 0.7 b = 0.7/0.2 = 3.5

    con lo que el polinomio queda:
    p(x) = 1 + 3.5 x - 2.5 x^2

    validando con python:
    tomando los puntos de prueba:

    xi = [ 0, 0.2, 0.4]
    fi = [ 1, 1.6, 2 ]
    >>>

    se obtiene la gráfica:

    se adjunta las instrucciones usadas para validar que el polinomio pase por los puntos requeridos.

    # 1Eva_IIT2017_T1 Aproximar a polinomio usando puntos
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    xi = [ 0, 0.2, 0.4]
    fi = [ 1, 1.6, 2 ]
    
    px = lambda x: 1 + 3.5*x - 2.5*(x**2)
    a = -0.5
    b = 1
    muestras = 21
    
    # PROCEDIMIENTO
    xj = np.linspace(a,b,muestras)
    pxj = px(xj)
    
    # SALIDA
    print(xj)
    print(pxj)
    
    # Gráfica
    plt.plot(xj,pxj,label='p(x)')
    plt.plot(xi,fi,'o', label='datos')
    
    plt.xlabel('x')
    plt.ylabel('f(x)')
    plt.grid()
    plt.legend()
    plt.show()
    
    

    Nota: Se puede intentar realizar los polinomios aumentando el grado, sin embargo cada término agrega un componente adicional a los términos anteriores por la forma (x - x0)k

    literal b

    se requiere el integral aproximado de f(x) en el intervalo limitado por los 3 puntos de la tabla:

    \int_{0}^{0.4}f(x) dx

    Esta aproximación con un polinomio es el concepto de integración numérica con la regla de Simpson de 1/3, tema desarrollado en la unidad 5

    I_{S13} = \frac{0.2}{3} \Big(1+4(1.6)+ 2 \Big) = 0.62666
  • s1Eva_IT2017_T4 Componentes eléctricos

    Ejercicio: 1Eva_IT2017_T4 Componentes eléctricos

    Desarrollo Analítico

    Solo puede usar las cantidades disponibles de material indicadas, por lo que las cantidades desconocidas de producción por componente se convierten en las incógnitas x0, x1, x2. Se usa el índice cero por compatibilidad con las instrucciones Python.

    Material 1 Material 2 Material 3
    Componente 1 5 x0 9 x0 3 x0
    Componente 2 7 x1 7 x1 16 x1
    Componente 3 9 x2 3 x2 4 x2
    Total 945 987 1049

    Se plantean las fórmulas a resolver:

    5 x0 +  7 x1 + 9 x2 = 945
    9 x0 +  7 x1 + 3 x2 = 987
    3 x0 + 16 x1 + 4 x2 = 1049
    

    Se reescriben en la forma matricial AX=B

    \begin{bmatrix} 5 & 7 & 9\\ 9 & 7 & 3 \\ 3& 16 & 4\end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix}= \begin{bmatrix} 945 \\ 987 \\ 1049 \end{bmatrix}

    Se reordena, pivoteando por filas para tener la matriz diagonalmente dominante:

    \begin{bmatrix} 9 & 7 & 3\\ 3 & 16 & 4 \\ 5& 7 & 9\end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix}= \begin{bmatrix} 987 \\ 1049 \\ 945 \end{bmatrix}

    Se determina el número de condición de la matriz, por facilidad con Python:

    numero de condicion: 4.396316324708121
    

    Obtenido con:

    # 1Eva_IT2017_T4 Componentes eléctricos
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    A = np.array([[9., 7.,3.],
                  [3.,16.,4.],
                  [5., 7.,9.]])
    
    B = np.array([987.,1049.,945.])
    # PROCEDIMIENTO
    
    # numero de condicion
    ncond = np.linalg.cond(A)
    
    # SALIDA
    print('numero de condicion:', ncond)
    

    Recordando que la matriz debe ser tipo real, se añade un punto a los dígitos.

    El número de condición es cercano a 1, por lo que el sistema si debería converger a la solución.

    Desarrollo con Python

    La forma AX=B permite usar los algoritmos desarrollados, obteniendo la solución. Se verifica el resultado al realizar la multiplicación de A con el vector respuesta, debe ser el vector B con un error menor al tolerado.

    respuesta con Jacobi
    [[62.99996585]
     [44.99998037]
     [34.9999594 ]]
    verificando:
    [[ 986.99943346]
     [1048.99942111]
     [ 944.99932646]]
    

    Si interpreta el resultado, se debe obtener solo la parte entera [63,45,35] pues las unidades producidas son números enteros.

    instrucciones:

    # 1Eva_IT2017_T4 Componentes eléctricos
    import numpy as np
    
    def jacobi(A,B,tolera,X0,iteramax=100, vertabla=False):
        ''' Método de Jacobi, tolerancia, vector inicial X0
            para mostrar iteraciones: tabla=1
        '''
        tamano = np.shape(A)
        n = tamano[0]
        m = tamano[1]
        diferencia = np.ones(n, dtype=float)
        errado = np.max(diferencia)
        X = np.copy(X0)
        xnuevo = np.copy(X0)
    
        itera = 0
        if vertabla==True:
            print('itera,[X0],errado')
            print(itera, xnuevo, errado)
        while not(errado<=tolera or itera>iteramax):
            
            for i in range(0,n,1):
                nuevo = B[i]
                for j in range(0,m,1):
                    if (i!=j): # excepto diagonal de A
                        nuevo = nuevo-A[i,j]*X[j]
                nuevo = nuevo/A[i,i]
                xnuevo[i] = nuevo
            diferencia = np.abs(xnuevo-X)
            errado = np.max(diferencia)
            X = np.copy(xnuevo)
            itera = itera + 1
            if vertabla==True:
                print(itera, xnuevo, errado)
        # Vector en columna
        X = np.transpose([X])
        # No converge
        if (itera>iteramax):
            X=itera
            print('iteramax superado, No converge')
        return(X)
    
    
    # INGRESO
    A = np.array([[9., 7.,3.],
                  [3.,16.,4.],
                  [5., 7.,9.]],dtype=float)
    
    B = np.array([987.,1049.,945.],dtype=float)
    tolera = 1e-4
    X0 = [0.,0.,0.]
    
    # PROCEDIMIENTO
    
    # numero de condicion
    ncond = np.linalg.cond(A)
    
    respuesta = jacobi(A,B,tolera,X0,vertabla=True)
    
    verifica = np.dot(A,respuesta)
    verificaErrado = np.max(np.abs(B-np.transpose(verifica)))
    
    # SALIDA
    print('numero de condicion:', ncond)
    print('respuesta con Jacobi')
    print(respuesta)
    print('verificar A.X:')
    print(verifica)
    print('max|A.X - B|')
    print(verificaErrado)
    

    al ejecutar el algoritmo se determina que se requieren 83 iteraciones para cumplir con con el valor de error tolerado.

  • s1Eva_IT2017_T2 Tanque esférico-volumen

    Ejercicio: 1Eva_IT2017_T2 Tanque esférico-volumen

    a. Planteamiento del problema

    V = \frac{\pi h^{2} (3r-h)}{3}

    Si r=1 m y V=0.75 m3,

    0.75 = \frac{\pi h^{2} (3(1)-h)}{3} 0.75 -\frac{\pi h^{2} (3(1)-h)}{3} = 0 f(h) = 0.75 -\frac{\pi h^{2} (3-h)}{3} = 0.75 -\frac{\pi}{3}(h^{2} (3)-h^{2}h) = 0.75 -\frac{\pi}{3}(3 h^{2} - h^{3}) f(h) = 0.75 -\pi h^{2} + \frac{\pi}{3}h^{3}

    b. Intervalo de búsqueda de raíz

    El tanque vacio tiene h=0 y completamente lleno h= 2r = 2(1) = 2, por lo que el intevalo tiene como extremos:

    [0,2]

    tanque esferico llenado altura h

    Verificando que exista cambio de signo en la función:

    f(0) = 0.75 -\pi (0)^{2} + \frac{\pi}{3}(0)^{3} = 0.75 f(2) = 0.75 -\pi (2)^{2} + \frac{\pi}{3}(2)^{3}= -3.4387

    y verificando al usar la gráfica dentro del intervalo:

    grafica tanque esferico llenado V vs h

    Tolerancia

    Se indica en el enunciado como 0.01 que es la medición mínima a observar con un flexómetro.

    tolera = 0.01

    cinta métrica


    c. Método de Newton-Raphson
    d. Método de Punto Fijo


    c. Método de Newton-Raphson

    El método de Newton-Raphson requiere la derivada de la función:

    x_{i+1} = x_i -\frac{f(x_i)}{f'(x_i)} f(h) = 0.75 -\pi h^{2} + \frac{\pi}{3}h^{3} f'(h) = -2\pi h + \pi h^{2}

    Tomando como punto inicial de búsqueda el extremo izquierdo del intervalo, genera una división para cero. Por lo que se mueve un poco a la derecha, algo más cercano a la raiz, viendo la gráfica por ejemplo 0.1

    x0 = 0.1

    iteración 1

    i =0

    f(0.1) = 0.75 -\pi (0.1)^{2} + \frac{\pi}{3}(0.1)^{3} =0.7196 f'(0.1) = -2\pi (0.1) + \pi (0.1)^{2} = -0.5969 x_{1} = x_0 -\frac{0.7496 }{-0.0625} = 1.3056 tramo = |x_{1}-x_{0}| = |0.1-1.3056 | = 1.2056

    iteración 2

    i =1

    f(1.3056) = 0.75 -\pi (1.3056)^{2} + \frac{\pi}{3}(1.3056)^{3} = -2.2746 f'(1.3056) = -2\pi (1.3056) + \pi (1.3056)^{2} =-2.8481 x_{1} = x_0 -\frac{-2.2746}{-2.8481} = 0.5069 tramo = |x_{2}-x_{1}|=|0.5069-1.3056|=0.7987

    iteración 3

    i =2

    f(0.5069) = 0.75 -\pi (0.5069)^{2} + \frac{\pi}{3}(0.5069)^{3} = 0.0789 f'(0.5069) = -2\pi (0.5069) + \pi (0.5069)^{2} =-2.3780 x_{1} = x_0 -\frac{0.0789}{-2.3780} = 0.5401 tramo = |x_{3}-x_{2}| =|0.5401 - 0.5069| = 0.0332

    Observe que el tramo disminuye en cada iteración , por lo que el método converge, si se siguen haciendo las operaciones se tiene que:

     [ xi, xnuevo, tramo]
    [[1.00000000e-01 1.30560920e+00 1.20560920e+00]
     [1.30560920e+00 5.06991599e-01 7.98617601e-01]
     [5.06991599e-01 5.40192334e-01 3.32007350e-02]
     [5.40192334e-01 5.39518667e-01 6.73667593e-04]]
    raiz 0.5395186666699257
    

    Instrucciones en Python

    para Método de Newton-Raphson

    # 1Eva_IT2017_T2 Tanque esférico-volumen
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    fx = lambda h: 0.75 - np.pi*(h**2)+(np.pi/3)*h**3
    dfx = lambda h: -2*np.pi*h+np.pi*(h**2)
    
    # Parámetros de método
    a = 0
    b = 2
    tolera = 0.01
    x0 = 0.1
    iteramax = 100
    
    # parametros de gráfica
    La = a
    Lb = b
    muestras = 21
    
    # PROCEDIMIENTO
    # Newton-Raphson
    tabla = []
    tramo = abs(2*tolera)
    xi = x0
    while (tramo>=tolera):
        xnuevo = xi - fx(xi)/dfx(xi)
        tramo = abs(xnuevo-xi)
        tabla.append([xi,xnuevo,tramo])
        xi = xnuevo
    tabla = np.array(tabla)
    
    # calcula para grafica
    hi = np.linspace(La,Lb,muestras)
    fi = fx(hi)
    gi = dfx(hi)
    
    # SALIDA
    print(' [ xi, xnuevo, tramo]')
    print(tabla)
    print('raiz', xnuevo)
    plt.plot(hi,fi)
    plt.plot(xi,fx(xi),'ro')
    plt.axhline(0, color='green')
    plt.xlabel('h')
    plt.ylabel('V')
    plt.title('Método Newton-Raphson')
    plt.show()
    

    Planteo con Punto Fijo


    d. Método de Punto Fijo

    Del planteamiento del problema en el literal a, se tiene que:

    0.75 = \frac{\pi h^{2} (3(1)-h)}{3}

    de donde se despeja una h:

    \frac{3(0.75)}{\pi (3(1)-h) } = h^{2} h = \sqrt{\frac{3*0.75}{\pi (3-h) }} h = \sqrt{\frac{2.25}{\pi (3-h) }}

    con lo que se obtienen las expresiones a usar en el método

    identidad = h g(h) = \sqrt{\frac{2.25}{\pi (3-h) }}

    El punto inicial de búsqueda debe encontrarse en el intervalo, se toma el mismo valor que x0 en el método de Newton-Raphson

    x0 = 0.10

    método punto fijo en tanque esférico

    Iteracion 1

    x_0= 0.10 g(0.10) = \sqrt{\frac{2.25}{\pi (3-(0.10) }}= 0.4969 tramo= |0.4969-0.10| = 0.3869

    Iteracion 2

    x_1= 0.4969 g(0.4969) = \sqrt{\frac{2.25}{\pi (3-(0.4969 ) }}= 0.5349 tramo= |0.5349- 0.4969| = 0.038

    Iteracion 3

    x_2 =0.5349 g(0.5349) = \sqrt{\frac{2.25}{\pi (3-(0.5349) }}= 0.5390 tramo= |0.5390 - 0.5349| = 0.0041

    con lo que se cumple el criterio de tolerancia, y se obtiene la raiz de:

    raiz = 0.5390

    Tabla de resultados, donde se observa que el tramo o error en cada iteración disminuye, por lo que el método converge.

     [i,xi,xi+1,tramo]
    [[1.         0.1        0.4969553  0.3969553 ]
     [2.         0.4969553  0.5349116  0.03795631]
     [3.         0.5349116  0.53901404 0.00410243]]
    raiz 0.5390140355891347
    >>> 
    

    Instrucciones en Python

    para Método de Punto-Fijo, recordamos que el método puede diverger, por lo que se añade el parámetro iteramax

    # 1Eva_IT2017_T2 Tanque esférico-volumen
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    fx = lambda h: h
    gx = lambda h: np.sqrt(2.25/(np.pi*(3-h)))
    
    a = 0.1
    b = 2
    tolera = 0.01
    iteramax = 100
    
    La = a
    Lb = b
    muestras = 21
    
    # PROCEDIMIENTO
    # Punto Fijo
    tabla = []
    i = 1 # iteración
    b = gx(a)
    tramo = abs(b-a)
    while not(tramo<=tolera or i>=iteramax):
        tabla.append([i,a,b,tramo])
        a = b
        b = gx(a)
        tramo = abs(b-a)
        i = i+1
    tabla.append([i,a,b,tramo])
    respuesta = b
    
    # Validar respuesta
    if (i>=iteramax):
        respuesta = np.nan
    tabla = np.array(tabla)
    
    # calcula para grafica
    hi = np.linspace(La,Lb,muestras)
    fi = fx(hi)
    gi = gx(hi)
    
    # SALIDA
    print(' [i, xi, xi+1, tramo]')
    print(tabla)
    print('raiz', respuesta)
    plt.plot(hi,fi, label = 'identidad', color='green')
    plt.plot(hi,gi, label = 'g(h)', color = 'orange')
    plt.plot(b,gx(b),'ro')
    plt.axhline(0, color='green')
    plt.xlabel('h')
    plt.ylabel('V')
    plt.title('Método Punto Fijo')
    plt.legend()
    plt.axhline(0, color='green')
    plt.show()
    

    Otra forma de probar la convergencia es que |g'(x)|<1 que se observa en la una gráfica adicional, lo que limita aún más el intervalo de búsqueda.

    Desarrollo en la siguiente clase.


  • s1Eva_IT2016_T3_MN Tasa interés anual

    Ejercicio: 1Eva_IT2016_T3_MN Tasa interés anual

    Propuesta de Solución  empieza con el planteamiento del problema, luego se desarrolla con el método de Bisección y método del Punto Fijo solo con el objetivo de comparar resultados.


    Planteamiento del problema

    La fórmula del enunciado para el problema es:

    A = P \frac{i(1+i)^{n}}{(1+i)^{n} -1}

    que con los datos dados se convierte a:

    5800 = 35000 \frac{i(1+i)^8}{(1+i)^8 -1} 35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800 =0

    que es la forma de f(x) = 0

    f(i)=35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800

    Intervalo de búsqueda

    Como el problema plantea la búsqueda de una tasa de interés, consideramos que:

    • Las tasas de interés no son negativas. ⌉(i<0)
    • Las tasas de interés no son cero en las instituciones bancarias (i≠0)
    • Las tasas muy grandes 1 = 100/100 = 100% tampoco tienen mucho sentido

    permite acotar la búsqueda a un intervalo (0,1].
    Sin embargo tasas demasiado altas tampoco se consideran en el problema pues el asunto es regulado (superintendencia de bancos), por lo que se podría intentar entre un 1% = 0.01l y la mitad del intervalo 50%= 0.5 quedando

    [0.01,0.5]

    Tolerancia

    si la tolerancia es de tan solo menor a un orden de magnitud que el valor buscado, se tiene que las tasas de interés se representan por dos dígitos después del punto decimal, por lo que la tolerancia debe ser menor a eso.

    Por ejemplo: tolerancia < 0.001 o aumentando la precisión

    tolera = 0.0001


    Método de la Bisección

    itera = 1

    a = 0.01, b = 0.5

    c = \frac{a+b}{2} = \frac{0.01+0.5}{2} = 0.255 f(0.01) = 35000 \frac{0.01(1+0.01)^8}{(1+0.01)^8-1}-5800 = -1225.83 f(0.5) = 35000 \frac{0.5(1+0.5i)^8}{(1+0.5)^8-1}-5800 = 12410.54

    con lo que se verifica que existe cambio de signo al evaluar f(x) en el intervalo y puede existir una raíz.

    f(0.255) = 35000 \frac{0.255(1+0.255)^8}{(1+0.255)^8-1}-5800 = 4856.70

    lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.255]

    El tramo permite estimar el error, reduce el intervalo a:

    tramo = b-a = 0.255-0.01 = 0.245

    valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

    itera = 2

    a = 0.01, b = 0.255

    c = \frac{0.01+0.255}{2} = 0.1325 f(0.01) = -1225.83 f(0.255) = 4856.70 f(0.1325 ) = 35000 \frac{0.1325(1+0.1325)^8}{(1+0.1325)^8-1}-5800 = 1556.06

    lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.1325]

    El tramo permite estimar el error, reduce el intervalo a:

    tramo = b-a = 0.1325-0.01 = 0.1225

    valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

    itera = 3

    a = 0.01, b = 0.1325

    c = \frac{0.01+0.1225}{2} = 0.071 f(0.01) = -1225.83 f(0.1325) = 1556.06 f(0.071 ) = 35000 \frac{0.071(1+0.071)^8}{(1+0.071)^8-1}-5800 = 89.79

    lo que muestra que f(x) tiene signos en a,c,b de (-) (+) (+), seleccionamos el intervalo izquierdo para continuar la búsqueda [0.01, 0.071]

    El tramo permite estimar el error, reduce el intervalo a:

    tramo = b-a = 0.071-0.01 = 0.061

    valor que todavía es más grande que la tolerancia de 10-4, por lo que hay que continuar las iteraciones.

    Para una evaluación del tema en forma escrita es suficiente para mostrar el objetivo de aprendizaje, el valor final se lo encuentra usando el algoritmo.

           raiz en:  0.06724243164062502
    error en tramo:  5.981445312500111e-05
    iteraciones:  13
    >>>

    Algoritmo en Python

    # 1Eva_IT2016_T3_MN Tasa interés anual
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    fx = lambda x: 35000*(x*(1+x)**8)/((1+x)**8 -1) -5800
    a = 0.01
    b = 0.5
    tolera = 0.0001
    
    # PROCEDIMIENTO
    cuenta = 0
    np.set_printoptions(precision=3)
    tramo = b-a
    while not(tramo<tolera):
        c = (a+b)/2
        fa = fx(a)
        fb = fx(b)
        fc = fx(c)
        cambia = np.sign(fa)*np.sign(fc)
        if cambia < 0: a = a b = c if cambia > 0:
            a = c
            b = b
        tramo = b-a
        cuenta = cuenta+1
    
    # SALIDA
    print('       raiz en: ', c)
    print('error en tramo: ', tramo)
    print('iteraciones: ',cuenta)
    


    Método del Punto Fijo

    El planteamiendo del punto fijo se realiza con x= g(x), por lo que se reordena la ecuación a:

    35000 \frac{i(1+i)^8}{(1+i)^8 -1}-5800 =0 35000 \frac{i(1+i)^8}{(1+i)^8 -1}=5800 \frac{i(1+i)^8}{(1+i)^8 -1}=\frac{5800}{35000} i=\frac{58}{350} \frac{(1+i)^8 -1} {i(1+i)^8}

    con lo que g(x) es:

    g(i)=\frac{58}{350} \frac{(1+i)^8 -1} {(1+i)^8}

    valor inicial de búsqueda

    Para el punto inicial i0 se puede usar uno de los extremos del intervalo propuesto en la sección de planteamiento. Para reducir aun más la búsqueda se pude seleccionar el punto intermedio

     i0 = 0.255

    Itera = 1

     i0 = 0.255

    g(0.255i)=\frac{58}{350} \frac{(1+0.255)^8 -1} {(1+0.255)^8} = 0.1387

    el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

    tramo = | nuevo-antes| = |0.1387 - 0.255| = 0.1163

    como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

    Itera = 2

    i1 = g(i0 ) = 0.1387

    g(0.1387)=\frac{58}{350} \frac{(1+0.1387)^8 -1} {(1+0.1387)^8} = 0.1071

    el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

    tramo = | nuevo-antes| = |0.1071 - 0.1387| = 0,0316

    como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

    Itera = 3

    i2 = g(i1 ) = 0.1071

    g(0.1071)=\frac{58}{350} \frac{(1+0.1071)^8 -1} {(1+0.1071)^8} = 0.0922

    el error se estrima como el tramo recorrido entre el valor nuevo y el valor inicial

    tramo = | nuevo-antes| = |0.0922 - 0.1071| = 0,0149

    como el tramo o error es aún mayor que el valor de tolera, se continúa con la siguiente iteración.

    Observación: Como el error disminuye entre cada iteración, se considera que el método converge, si se realizan suficientes iteraciones se cumplierá con el valor de tolerancia y se habrá llegado a la precisión requerida.

    Usando el algoritmo se tiene:

    iteraciones: 17
        raiz en: 0.06754389199556779
    >>>

    Algoritmo en Python

    # Algoritmo de Punto Fijo
    # x0 valor inicial de búsqueda
    # error = tolera
    
    import numpy as np
    def puntofijo(gx,antes,tolera, iteramax=50):
        itera = 1 # iteración
        nuevo = gx(antes)
        tramo = abs(nuevo-antes)
        while(tramo>=tolera and itera<=iteramax ):
            antes = nuevo
            nuevo = gx(antes)
            tramo = abs(nuevo-antes)
            itera = itera+1
        respuesta = nuevo
        # Validar respuesta
        if (itera>=iteramax ):
            respuesta = np.nan
        print('iteraciones:',itera)
        return(respuesta)
    
    # PROGRAMA ---------
    
    # INGRESO
    gx = lambda i: (58/350)*((1+i)**8-1)/((1+i)**8)
    
    a = 0       # intervalo
    b = 0.5
    
    x0 = 0.255
    tolera = 0.0001
    iteramax  = 50      # itera máximo
    
    # PROCEDIMIENTO
    respuesta = puntofijo(gx,0.255,tolera,iteramax)
    
    # SALIDA
    print('    raiz en:',respuesta)
    
  • s1Eva_IT2015_T4 Lingotes metales

    Ejercicio: 1Eva_IT2015_T4 Lingotes metales

    Se plantea que cada lingote debe aportar una proporción xi al lingote nuevo a ser fundido.

    Se dispone de los compuestos de cada lingote por filas:

    compuesto = np.array([[ 20, 50, 20, 10],
                          [ 30, 40, 10, 20],
                          [ 20, 40, 10, 30],
                          [ 50, 20, 20, 10]])
    proporcion = np.array([ 27, 39.5, 14, 19.5])

    El contenido final de cada componente basado en los aportes xi de cada lingote para cada componente.

    Ejemplo para los 27 gramos de oro

    20x_1 + 30x_2+ 20x_3 + 50x_4 = 27

    se realiza el mismo procedimiento para los otros tipos de metal.

    50x_1 + 40x_2+ 40x_3 + 20x_4 = 39.5 20x_1 + 10x_2+ 10x_3 + 20x_4 = 14 10x_1 + 20x_2+ 30x_3 + 10x_4 = 19.5

    Las ecuaciones se escriben en la forma matricial Ax=B

    \begin{bmatrix} 20 && 30&& 20 &&50 \\ 50 && 40 && 40 && 20 \\ 20 && 10 && 10 && 20 \\ 10 && 20 && 30 && 10 \end{bmatrix} = \begin{bmatrix} x_1 \\ x_2 \\x_3 \\ x_4 \end{bmatrix} = \begin{bmatrix} 27 \\ 39.5 \\ 14 \\ 19.5 \end{bmatrix}

    Para resolver se plantea la matriz aumentada

    \begin{bmatrix} 20 && 30&& 20 &&50 && 27\\ 50 && 40 && 40 && 20 && 39.5\\ 20 && 10 && 10 && 20 && 14 \\ 10 && 20 && 30 && 10 && 19.5 \end{bmatrix}

    se pivotea por filas la matriz:

    \begin{bmatrix} 50 && 40 && 40 && 20 && 39.5\\ 20 && 30&& 20 &&50 && 27\\ 10 && 20 && 30 && 10 && 19.5 \\ 20 && 10 && 10 && 20 && 14 \end{bmatrix}

    Para eliminar hacia adelante:

    \begin{bmatrix} 50 && 40 && 40 && 20 && 39.5 \\ 20 - 50\frac{20}{50} && 30-40\frac{20}{50} && 20-40\frac{20}{50} && 50-20\frac{20}{50} && 27-39.5\frac{20}{50}\\ 10 - 50\frac{10}{50} && 20-40\frac{10}{50} && 30-40\frac{10}{50} && 10-20\frac{10}{50} && 19.5-39.5\frac{10}{50} \\ 20 - 50\frac{20}{50} && 10-40\frac{20}{50} && 10-40\frac{20}{50} && 20-20\frac{20}{50} && 14-39.5\frac{20}{50} \end{bmatrix}

    continuando con el desarrollo:

    Elimina hacia adelante
    [[50.  40.  40.  20.  39.5]
     [ 0.  14.   4.  42.  11.2]
     [ 0.  12.  22.   6.  11.6]
     [ 0.  -6.  -6.  12.  -1.8]]
    Elimina hacia adelante
    [[ 50.  40.  40.      20.  39.5]
     [  0.  14.   4.      42.  11.2]
     [  0.   0.  18.5714 -30.   2. ]
     [  0.   0.  -4.2857  30.   3. ]]
    Elimina hacia adelante
    [[ 50.  40.  40.      20.     39.5   ]
     [  0.  14.   4.      42.     11.2   ]
     [  0.   0.  18.5714 -30.      2.    ]
     [  0.   0.   0.      23.0769  3.4615]]
    Elimina hacia adelante
    [[ 50.  40.   40.      20.     39.5  ]
     [  0.  14.    4.      42.     11.2  ]
     [  0.   0.   18.5714 -30.      2.   ]
     [  0.   0.    0.      23.0769  3.4615]]
    Elimina hacia atras
    [[ 1.    0.    0.    0.    0.25]
     [ 0.    1.    0.    0.    0.25]
     [ 0.    0.    1.   -0.    0.35]
     [ 0.    0.    0.    1.    0.15]]
    el vector solución X es:
    [[0.25]
     [0.25]
     [0.35]
     [0.15]]
    verificar que A.X = B
    [[39.5]
     [27. ]
     [19.5]
     [14. ]]
    

    Las proporciones de cada lingote a usar para el nuevo lingote que cumple con lo solicitado son:

    [0.25, 0.25, 0.35, 0.15]


    Algoritmo en python

    usado para la solución es:

    # Método de Gauss-Jordan ,
    # Recibe las matrices A,B
    # presenta solucion X que cumple: A.X=B
    import numpy as np
    import matplotlib.pyplot as plt
    
    # INGRESO
    A = np.array([[ 50., 40, 40, 20],
                  [ 20., 30, 20, 50],
                  [ 10., 20, 30, 10],
                  [ 20., 10, 10, 20]])
    B1 = np.array([ 39.5, 27, 19.5, 14])
    
    B = np.transpose([B1])
    
    # PROCEDIMIENTO
    casicero = 1e-15 # 0
    AB = np.concatenate((A,B),axis=1)
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    
    print('matriz aumentada: ')
    print(AB)
    # Gauss elimina hacia adelante
    # tarea: verificar términos cero
    for i in range(0,n,1):
        pivote = AB[i,i]
        adelante = i+1 
        for k in range(adelante,n,1):
            if (np.abs(pivote)>=casicero):
                factor = AB[k,i]/pivote
                AB[k,:] = AB[k,:] - factor*AB[i,:]
        print('Elimina hacia adelante')
        print(AB)
    
    # Gauss-Jordan elimina hacia atras
    ultfila = n-1
    ultcolumna = m-1
    for i in range(ultfila,0-1,-1):
        # Normaliza a 1 elemento diagonal
        AB[i,:] = AB[i,:]/AB[i,i]
        pivote = AB[i,i] # uno
        # arriba de la fila i
        atras = i-1 
        for k in range(atras,0-1,-1):
            if (np.abs(AB[k,i])>=casicero):
                factor = pivote/AB[k,i]
                AB[k,:] = AB[k,:]*factor - AB[i,:]
            else:
                factor= 'division para cero'
    print('Elimina hacia atras')
    print(AB)
    
    X = AB[:,ultcolumna]
    
    # Verifica resultado
    verifica = np.dot(A,X)
    
    # SALIDA
    print('el vector solución X es:')
    print(np.transpose([X]))
    
    print('verificar que A.X = B')
    print(np.transpose([verifica]))
    

    Tarea: Revisar sobre la última pregunta.

  • s1Eva_IT2015_T2 Salida cardiaca

    Ejercicio: 1Eva_IT2015_T2 Salida cardiaca

    Solución presentada como introducción al tema de interpolación y solución de sistemas de ecuaciones.
    No realiza el literal c, no se desarrolla el tema de integrales.

    Note que el desarrollo del tema permite aumentar el grado del polinomio de interpolación, lo que se afecta al tamaño del sistema de ecuaciones (matriz).

    Los valores obtenidos con la solución propuesta son:

    solución para X por Gauss-Seidel
    [[-0.42175867]
     [ 0.15610893]
     [ 0.02736763]]
    verifica que A.X = B:
    [[ -7.02831455e-05]
     [  1.50012970e+00]
     [  3.20000000e+00]]
    polinomio interpolación, con puntos:  3
    0.0273676337623498*x**2 + 0.156108926206362*x - 0.421758670607596
    >>>
    

    La gráfica para observar los datos experimentales es:
    s1EIT2015 salida cardiaca 01

    La gráfica con polinomio de interpolación de grado 2, con tres puntos:

    s1EIT2015 salida cardiaca 02

    instrucciones del problema para la solución por partes en python:

    # 1ra Evaluación I Término 2015
    # Tema 2. Flujo de sangre en corazón
    # Tarea: parte c), no se ha realizado el áre bajo curva
    #        falta calcular salida cardiaca.
    import numpy as np
    import matplotlib.pyplot as plt
    
    # Gráfica de datos experimentales:
    t = np.array([2,6,9,12,15,18])
    y = np.array([0,1.5,3.2,4.1,3.4,2.0])
    
    # SALIDA
    plt.plot(t,y)
    plt.title('datos del experimento: t vs concentración ')
    plt.show()
    
    # Sistema de ecuaciones para aproximar a polinomio grado 2
    # para grado dos usa tres puntos,
    # por ejemplo usando el punto [ 2, 0] del experimento
    # a + b*(2) + c*(2**2) =  0
    A = np.array([[1,2,2**2],
                  [1,6,6**2],
                  [1,9,9**2]])
    B = np.array([[0],
                  [1.5],
                  [3.2]])
    
    tolera = 0.0001
    X = np.zeros(len(B), dtype = float)
    
    # usando numpy para solucion de matrices
    # Xnp = np.linalg.solve(A,B)
    # print('solución para A.X=B con numpy')
    # print(Xnp)
    
    # algoritmo Gauss-Seidel
    iteramax=100
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    
    itera = 0
    while (errado>tolera or itera>iteramax):
        for i in range(0,n,1):
            nuevo = B[i]
            for j in range(0,m,1):
                if (i!=j): # excepto diagonal de A
                    nuevo = nuevo-A[i,j]*X[j]
            nuevo = nuevo/A[i,i]
            diferencia[i] = np.abs(nuevo-X[i])
            X[i] = nuevo
        errado = np.max(diferencia)
        itera = itera + 1
    # Vector en columna
    X =  np.transpose([X])
    # No converge
    if (itera>iteramax):
        X=0
    
    Xgs = X
    
    # Metodo numérico Gauss_Seidel
    verifica = np.dot(A,Xgs)
    print('solución para X por Gauss-Seidel')
    print(Xgs)
    
    # verificar resultado
    print('verifica que A.X = B: ')
    print(verifica)
    
    # Observar interpolacion con polinomio creado
    pt = lambda t: Xgs[0,0]+ Xgs[1,0]*t + Xgs[2,0]*t**2
    
    ti = np.linspace(2,18,501)
    pti = pt(ti)
    
    plt.plot(ti,pti, label = 'interpolacion')
    plt.plot(t,y,'*', label = 'datos experimento')
    plt.title('interpolación con polinomio')
    plt.legend()
    plt.show()
    
    # polinomio en sympy
    import sympy as sp
    x = sp.Symbol('x')
    polinomio = 0
    k = len(Xgs)
    for i in range(0,k,1):
        polinomio = polinomio + Xgs[i,0]*x**i
    print('polinomio interpolación, con puntos: ', k) 
    print(polinomio)